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The molecular component of the Galaxy is comprised of turbulent, 
magnetized clouds, many of which are self-gravitating and form stars. To 
develop an understanding of how these clouds' kinetic and structural evolution 
may depend on their level of turbulence, mean magnetization, and degree of 



O 



self-gravity, we perform a survey of direct numerical MHD simulations in which 
three parameters are independently varied. Our simulations consist of solutions 
to the time- dependent MHD equations on a two-dimensional grid with periodic 
q ■ boundary conditions; an additional "half" dimension is also incorporated 

+3 • as dependent variables in the third Cartesian direction. Two of our survey 

parameters - the mean magnetization parameter f3 = (? smm J^\it ven and the 
Jeans number rtj = L c i ou d/ Lj eans - allow us to model clouds which either meet 
^ . or fail conditions for magneto- Jeans stability and magnetic criticality. Our third 

survey parameter - the sonic Mach number M. = <J V eiocit y / c soun d - allows us 
to initiate turbulence of either sub- or super- Alfvenic amplitude; we employ 
an isothermal equation of state throughout. We evaluate the times for each 
cloud model to become gravitationally bound, and measure each model's kinetic 
energy loss over the fluid flow crossing time. We compare the evolution of 
density and magnetic field structural morphology, and quantify the differences in 
the density contrast generated by internal stresses, for models of differing mean 
magnetization. We find that the values of (3 and rij, but not the initial Mach 
number A4, determine the time for cloud gravitational binding and collapse: 
for mean cloud density n# 2 = 100cm -3 , unmagnetized models collapse after 
~ 5Myr, magnetically supercritical models generally collapse after 5 — lOMyr 
(although the smallest magneto- Jeans stable clouds survive gravitational collapse 
until t ~ 15Myr), while magnetically subcritical clouds remain uncollapsed over 
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the entire simulations; these cloud collapse times scale with the mean density as 

— 1/2 

t g oc n H2 . We find, contrary to some previous expectations, less than a factor of 
two difference between turbulent decay times for models with varying magnetic 
field strength; the maximum decay time, for B ~ 14/iG and n# 2 = 100cm -3 , is 
1.4 flow crossing times t cross = L/a ve i ocit y (or 8Myr for typical GMC parameters). 
In all models, we find turbulent amplification in the magnetic field strength up 
to at least the level f3 pert = c 2 soun dl & v \ifven = 0-1) with the turbulent magnetic 
energy between 25-60% of the turbulent kinetic energy after one flow crossing 
time. We find that for non-self-gravitating stages of evolution, when clouds have 
M. = 5 — 10, the mass-averaged density contrast magnitudes (log(p/p)) are in 
the range 0.2 — 0.5, with the contrast increasing both toward low- and high- 
(3. Although our conclusions about density statistics may be affected by our 
isothermal assumption, we note that only the more strongly-magnetized models 
appear consistent with estimates of clump/interclump density contrasts inferred 
in Galactic GMCs. 



Introduction 



All present-day star formation in our Galaxy is observed to take place in cold molecular 
clouds. These clouds are highly turbulent (with sonic Mach numbers 10 or more for the most 
massive clouds, e.g. piitz (1993)| ), and appear to be quite inhomogeneous based on both 
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, Lada, Alves, & 





Lada ( 1998 )|) . Because magnetic field measurements are difficult to obtain, the distributions 
of magnetic field strengths and directions within clouds are less well characterized than 
other cloud properties; the mean field strengths are probably no less than a few /iG, and no 
more than a few tens of /iG (e.g. [Heiles et al. (1993 )| , |Grutcher et al (1993 )| , [Ireland et al| 
(1996)| , |Urutcher (1998)| ), and the spatial power spectrum of the magnetic field is probably 



dominated by components with close to the largest scales possible, within structures of a 



given density regime (e.g. |Goodman fc Heiles ( 199471 , |Goodman et al (1995)| , |Schleuning 



1998)|) Given the extremely low thermal pressure in molecular clouds, magnetic stresses 



are expected to be at least as important as gas pressure in governing cloud evolution, and 
may be much more so. In particular, since sufficiently strong mean magnetic fields would 
render the internal motions sub-Alfvenic and the cloud self-gravity subcritical, it has long 
been argued that magnetization may significantly inhibit turbulent dissipation and prevent 
gravitational collapse ( |Arons fc Max (1975)| , |Mouschovias fc Spitzer (1976)| , |Shu et aT. 
P7| [CTcKee et al. (1993)|). 
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The nonlinear amplitudes and structural irregularity of the turbulence seen in molecular 
clouds require time-dependent numerical simulations in order to model their evolution 
theoretically. This is the second in a series of papers reporting on the results of numerical 
experiments designed to characterize how magnetic fields at a range of strengths affect the 
structure and evolution of cold, turbulent, self-gravitating clouds. In Paper I QGammie fc 



Ostriker (1996)|) , we outlined (see also |Ostriker ( 1997)1 ) the observational situation, the 



theoretical background on the nature of ideal turbulent flows and applications to cold 
Galactic clouds, and the issues that can be addressed by numerical simulations; we then 
presented results of the first extensive survey (in "1 2/3D" slab geometry, i.e. with one 
independent spatial variable, but all three components of dependent spatial variables) of 
model cloud evolution in the ideal MHD limit. The survey considered clouds both with 
and without gravity, and either with an input spectrum of turbulence that decayed over 
time, or with an ongoing excitation of turbulence through velocity perturbations introduced 
at spatial scales smaller than the simulation box. Among the conclusions of Paper I was 
the finding that, when sufficient initial energy in the form of transverse disturbances is 
input to a cloud model and then evolved, gravitational collapse along the mean field can 
be suppressed for times > t g (see eq. [| below). Paper I also reported on the large density 
contrasts, and correspondingly small filling factors, which arise in freely-decaying turbulent 
MHD flows (e.g. a typical model having 50% of the mass at density larger than three 
times the mean, and 10% of the mass at nine times greater than the mean). Because these 
conclusions could potentially be affected by the restricted geometry of the slab-symmetric 
models - where no spatial variations perpendicular to the mean field are possible - we 
re-examine these issues in the present paper with a less restricted geometry. 

Recent numerical work by other groups has addressed a variety of related questions. 
Elmegreen (1997")| has found that hierarchical density structure arises in ID MHD cloud 



models which are subject to mechanical forcing at their boundaries. [Passot, Vazquez- 



Semadeni, fc Pouquet (1995)| and |Vazquez-Semadeni, Passot, fc Pouquet (1996)| have 



modeled cloud formation and evolution in the ISM on ~ kpc scales, incorporating a number 
of physical effects (Galactic rotation and shear, heating from localized star formation, 
and cooling) in addition to magnetic fields. Padoan fc Nordlund (1997] and MacLow et 



|al ( 1998^1 have simulated the decay of Mach-5 turbulence in 3D cloud models with weak 
(c 2 s /v\ = 2, 1, respectively) and moderate (c 2 s /v\ = 0.02,0.04, respectively) magnetic fields, 
and find less than a factor 2 difference in the kinetic energy decay for the two cases. 

Because the matter distributions in clouds should eventually influence star formation 
and the IMF, numerical modelers have taken a particular interest in analyzing the density 
structure of their simulated clouds. In particular, there has been considerable recent 
attention devoted to the shape of the density distribution function obtained in simulations 
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of compressible turbulence. Based on a variety of non-self-gravitating simulations at varied 
Mach numbers, with/without turbulent forcing, and using a variety of polytropic indices 
7, two different groups flPadoan, Jones, fc INordlund (1997)| , |JNordlund fc Padoan (1998)| ; 



Vazquez-Semadeni (1994)| , [Scalo et al (1998)] , and [Passot fc Vazquez-Semadeni (1998]]) have 



concurred in finding that log- normal density distributions are expected when 7 = 1, and 
the latter group finds skewed log- normal/power law distributions are expected when 7^1. 
These groups differ considerably, however, in their predictions of the scaling of the mean 
density contrast magnitude with Mach number when 7 = 1 (see §3.3). 

In this paper, we extend our previous slab-geometry survey (Paper I) by performing 
a survey in 2 1/2 dimensions^ of the evolution of turbulence in model clouds with a 
range of sizes, Mach numbers, and magnetizations. In particular, since incompressible 
motions (i.e. nonlinear Alfven waves) have often been considered as a weakly-dissipative 
effective "turbulent" pressure that can prevent cloud collapse over several free-fall times 
(see Paper I), we focus on the evolution of self-gravitating model clouds in which an initial, 
incompressible turbulent spectrum is initiated, but no subsequent small-scale energy is 
injected. Thus, although our models are highly idealized and should not be thought of 
as representing real clouds, our results are a useful guide for understanding how cloud 
dynamics would proceed in the absence of potential turbulent inputs, many of which have 
been attributed to the onset of internal star formation (e.g. [Norman fc Silk (1980) , Micsch 



fc Bally (1994)| , [Passot, Vazquez-Semadeni, fc Pouquet (1995)| ). The plan of the paper 



is as follows: §2 describes our numerical experiments, including the computational code 
we use (§2.1), the definitions of relevant dimensionless parameters (§2.2), and our choices 
of parameters for the survey (§2.3). In §3, we present our results on the evolution of 
cloud energies (§3.1), structural appearance (§3.2), and density distributions (§3.3) in our 
models. We summarize our results and discuss the principal implications of our work for 
the evolution and structure of observed clouds in §4. 



2. Simulations 

2.1. Numerical Method 

To create the model cloud evolutions presented in this paper, we integrate the 
equations of ideal, compressible, self-gravitating MHD using a variant of the ZEUS 



1 i.e. two spatial independent variable, but three spatial components to all dependent variables; this 
allows spatial variations both parallel and perpendicular to the mean magnetic field. 



- 5 - 



algorithm developed by [Stone fc Norman (1992a) , Stone fc Norman (1992b) . ZEUS is an 



operator-split, finite difference method on a staggered mesh. The magnetic field is evolved 
using "constrained transport" ( [Evans fc Hawley (1988)[ ), which guarantees that V • B = to 



machine precision. The transverse components of the magnetic field are evolved using the 
"method of characteristics" which guarantees accurate propagation of Alfvenic disturbances. 
The ZEUS algorithm has been extensively tested and used in a wide variety of astrophysical 
applications. In the present models, dissipation occurs on small scales in part due to purely 
numerical effects associated with discretizing the MHD equations (i.e. from zone-to-zone 
averaging of velocity, magnetic field, or density variations), and in part due to an artificial 
viscosity introduced to capture shocks. We do not include ambipolar diffusion (which is 
important on small scales in GMCs), or other explicit resistivity or viscosity. 

Our implementation of ZEUS is "2 1/2 dimensional" (2.5 D). This peculiar 
nomenclature means that there are two independent spatial variables (x,y), but that all 
(three) components of the velocity and magnetic field vectors are still evolved as dependent 
variables. In this sense, there are two fully dynamical dimensions and a third dimension in 
which the flow is symmetric. The boundary conditions in all models are periodic. 

The gravitational potential is determined from V 2 = 4irG(p — p); the mean density 
appears because of the periodic boundary conditions. We solve for the gravitational 
potential using the Fourier method. The gravitational kernel is 



k = 2vrGp k ((1 - cos(k x Ax))/Ax 2 + (1 - cos(k y Ay))/Ay 2 ) . (1) 

This ensures that the discrete representation of Poisson's equation is satisfied. 

We adopt, for all models, an isothermal equation of state, meaning that p = c^p, where 
c 2 is a constant in both space and time. Thus we do not integrate an internal energy 
equation. A non-adiabatic assumption is a physically reasonable approximation for the 
higher- density (ra# 2 ~ 10 2 — 10 4 cm~ 3 ) portions within molecular clouds because the cooling 
time is very short, leading to a range of observed temperatures 10 — 20K over molecular 
regions not directly associated with high-luminosity star formation (e.g. goldsmith (1987)] 
and references therein). If radiative processes are sufficiently rapid compared to dynamical 
times, then an effective polytropic index 7 may be found by balancing heating and cooling 
rates locally ( [Scalo et al (1998]]) . The value of the effective 7 increases from 7 < 1 to 
7 ^ 1 near nn 2 ~ 10 3 cm -3 , with the increase due to the inefficiency of cooling in optically 
thick regions (radiative trapping). The details of the transition in this effective 7 depend 
on the relation of the density at a given velocity to the column density at that velocity 
(cf. goldsmith fc hanger (1978]| , |de Jong, Dalgarno, fc Boland (1980)1 |Neufeld, Lep"p71i| 



|Mclnick (1995)| ); this relation however remains statistically uncertain (as well as spatially 
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multivalued) in a turbulent medium. Since no single choice of an effective 7 applies to 
all density regimes of interest in the interior of a GMC, it would be highly desirable to 
perform numerical calculations in which both thermodynamics and MHD equations are 
evolved in detail. This would be a formidable calculation to perform realistically, because 
the molecular cooling lines are optically thick; to date, no group has done hydrodynamic 
calculations with even model heating and cooling included, in the regime appropriate for 
GMCs. At this stage, when the effects of magnetic fields on turbulent dynamics are not 
understood, we view the choice 7 = 1 as an acceptable compromise because most of the gas 
evolves to higher-than-average densities for which the temperature likely varies by less than 
a factor ~ 2, even if a strictly isothermal effective 7 applies only to n H2 ~ 10 3 cm -3 . The 
isothermal assumption is less accurate for the low-density "interclump" regions that form 
as a result of turbulent dynamics. We comment further on the implications of our choice of 
7 in §3.3. 

For all the simulations, we choose the initial velocity field to be locally divergence-free 
(i.e. the individual Fourier components are constrained to satisfy Vk ■ k = 0), and isotropic 
in direction. The velocity amplitude is drawn from a Gaussian distribution with power 
spectrum (|vk| 2 ) oc A; -3 , for k = \k x x + k y y\ and 27r/L < k{ < 128(2ir/L). The initial 
kinetic energy therefore obeys E = T,^^ where E^ = (l/2)v^ oc k~ 3 , corresponding to a 
relationship <r v (R) oc R 1 ^ 2 between size scale R and the velocity dispersion a v averaged over 
that scale; we choose this initial spectrum to match the linewidth-size relation for observed 
self-gravitating, molecular clouds on a variety of scales. For comparison, we note that a 
Kolmogorov spectrum in 2D has (|v k | 2 ) oc k~ 8 ^ 3 , so that a v (R) oc _R 1//3 . 

Initially, the density for all models is taken to be uniform; all inhomogeneities are 
generated self-consistently from divergences in the velocity field created by internal pressure 
forces, gravity, and magnetic stresses (since V • v = initially). The initial magnetic field 
is also taken to be uniform, pointing in the x direction (right-left in the model snapshots 
shown). Notice that the uniformity of B and p in the initial conditions implies that the 
mass-to-flux ratio along each field line is constant. In snapshots shown, the symmetric 
direction z points out of page. 

We have tested our implementation of the ZEUS algorithm using a variety of test 
problems: linear waves, standard advection tests, shock tubes and the parametric instability 
of circular polarized Alfven waves. Of particular interest for what follows is whether the 
simulations are converged at our standard resolution, 256 2 zones. We have tested this 
by evolving an initial Gaussian random velocity field, as in all our decay experiments, 
while holding the amplitude and phase of the field's Fourier components fixed as the 
resolution is increased from 32 2 to 512 2 zones. The convergence test models use n.j = 2.5, 



- 7- 



EK,init = 50, and (3 = 0.02 (see §2.2 for definitions). We find a secular decrease in the ratios 
E 2n /E n of turbulent energies (= perturbed magnetic and kinetic energy) at successive 
resolutions; at time 0.1t s these ratios decrease from E GA / E 32 = 1.19 to E 128 /E M = 1.11 to 
E 256 /E 128 = 1.08 to E 512 /E 256 = 1.06. We also find only a 5% difference in the gravitational 
binding time (see §3.1) between the 256 2 and 512 2 models. Thus, at a resolution of 256 2 
zones, the evolution is satisfactorily converged. 



2.2. Scalings and definitions 

The fundamental dimensional units for our simulations are the length L of the box 
edge, the mean density p = Mj L 3 of matter in the box, and the (isothermal) sound speed 
c s = yJkT/fj,. Because self-gravity is crucial in the problems under study, it is useful to 
describe the size of the simulated region relative to the minimum Jeans-unstable wavelength 
Lj = c s (7r /Gp) 1 / 2 at the cloud's mean density and temperature. For each simulation we 
therefore define the size via the Jeans number nj = L/Lj; using c s = 0.19kms _1 (T/10K) 1 / 2 
for // = 2Am p the corresponding physical size is 

L = njLj = nj x 1.9pc ( —) ' " ( ^ , ) . (2) 
J J J y U0K7 V10 2 cm- 3 / v ; 

With the parameter values nj = 2, 3, 4 used in the present simulations, the box edge would 
respectively correspond to approximately 4, 6, or 8pc at the mean conditions in a molecular 
cloud {tih 2 ~ 100), a factor 3 smaller for conditions within a cloud clump (n# 2 ~ 10 3 ), or a 
factor 10 smaller for conditions in a dense core (n# 2 ~ 10 4 ). 

Time in our numerical experiments is measured either in units of the sound-crossing 
time t s = L/c s for the simulated region, or in terms of a characteristic gravitational 
contraction time t g = Lj/c s = (n/Gp) 1 / 2 !^ The corresponding physical time depends on the 
mean density according to 

% 2 ^ 2 



= = 9.9 Myr(_^j . (3) 

Thus for rij = 2,3, or 4, the interval 0.1t s would respectively correspond to approximately 
2, 3, or 4 Myr at the typical mean density in a molecular cloud, and a factor 3 or 10 shorter 
at clump or core densities; the interval 0.1t g ~ lMyr at typical mean GMC conditions. 



2 For reference, the free-fall collapse times for a cold sheet, cylinder, or sphere are, respectively 0.25t g , 
0.28t gi and 0.3lt g . 
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The total mass within a cube with faces the same size as the simulated region 
corresponds to M — n 3 jMj, where the Jeans mass at the mean conditions is 

Mj = pL] = 48 M ( ^ X 1 ' 2 ( -^-) 3/2 . (4) 

J H J ° vio 2 cm- 3 / viok; v ; 

At the beginning of all the simulations, the magnetic field is uniform with strength Bq 
and points in the x direction. The level of the mean magnetic field in each simulation is 
characterized by the ratio 



~2 J2 



P ~ v\ ~ B 2 /(4np) ~ °-° 21 (lOK) (lOW^) (lO/lc) (5) 

1 /p / rp \ 1/2 / n „ \ 1/2 

corresponding to B = lApG(3~ 1 ( Iok j ( lo^cm- 3 J ' decreasing f3 corresponds 
to increasing the importance of magnetic fields to the dynamics. For the values 
j3 = 0.01, 0.1, 1.0 used in the present set of simulations, the mean magnetic field is 
B = 14, 4.4, and 1.4/xG respectively, for fiducial values T = 10K and n# 2 = 100cm -3 . 

Uniform clouds are unstable to compressions transverse to the mean magnetic 
field ("magneto- Jeans unstable") when the magnetosonic wave crossing time exceeds t g 
flChandrasekhar fc Fermi (1953)[) ; for v\ >> c s this occurs when nj > /? _1//2 , corresponding 
to having the cross-field column density Nh 2 > 4.1 x 10 21 cm~ 2 (I?o/10/iG) . A closely- 
related criterion is used to evaluate a more stringent stability condition, as follows: If a 
uniform cubic cloud with size L = njLj condenses along the mean field to make a cold 
"pancake" (with central density p = E 2 7rG/(2c 2 ) = (7mj) 2 p/2), the pancake will have 
ratio of surface density to magnetic field S/5 = (rij/3 1 / 2 ) / (2G 1 ^ 2 ). This magnetized 
pancake is unstable to fragmentation when Tj/Bq > l/(27rG 1 / 2 ) (corresponding to 
Nh 2 > 1.3 x 10 21 cm _2 (i?o/10/iG)), so that the corresponding instability criterion in terms 
of the original cloud's parameters is nj > /5 _1 / 2 7r _1 . Clouds with nj > (3~ l l 2 ix~ l are termed 
"supercritical;" clouds with the inequality reversed are termed "subcritical." We note 
that the gravitational plus magnetic energy of the "pancake" is still negative provided 
6/(7mj) 4 < f3. 

In terms of the fundamental units, all energies in the problem are given in units of 
pc 2 L 2 ; because the z direction is symmetric, this is the energy per unit length dz. For a 
cubic cloud, the corresponding total energy unit would be pc 2 L 3 . Where the kinetic energy 



-2 



3 The same instability criterion applies when considering either the ratio T,/Bq for a uniform magnetized 
sheet or the central mass-to- magnetic- flux ratio of a magnetized cloud (Tomisaka, Ikeuchi, & Nakamura 
(1988)|). 
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is given, it corresponds to 

1 rL/2 rL/2 

E K = - / dxdy(vl + v 2 + v 2 z )p; (6) 

2 J-L/2 J-L/2 y 

the square of the Mach number M. = cr v /c s (where a v represents the total velocity 
dispersion) is therefore J\A 2 = (<j v /c s ) 2 = 2[Ek/ (pL 2 c 2 )}, i.e. twice the normalized kinetic 
energy. The total magnetic energy corresponds to 



1 r L/2 r L/2 

Eb = - I / dxdy(B 2 x + B 2 + B 2 ); (7) 

07T J-L/2 J-L/2 y 



for the energy in the perturbed magnetic field, SEb, we subtract out the mean- field energy 
1/(2/?) from equation (|7|). Because we assume periodic boundary conditions, the mean 
magnetic field (corresponding to the k = Fourier component of B), and hence the 
mean-field energy, remains unchanged over the course of any simulation. The gravitational 
energy is 

1 rL/2 f L/2 

E G = - / dxdy{(j>p) (8) 

2 J-L/2 J-L/2 

where the gravitational potential is computed from the modified Poisson equation in periodic 
boundary conditions (§1.1). We note that because equation ([[]) implies (fi/c 2 oc L 2 Gp/c 2 , 
the scaled gravitational energy will obey Eg/ (pL 2 c 2 ) oc n 2 . 



2.3. Survey Model Parameters 

As described in §1, observed clouds present a range of properties, with some more 
tightly constrained than others. From a theoretical point of view, the observed bulk cloud 
properties of size, mean magnetization, and velocity dispersion can be described in terms 
of the dimensionless parameters nj, j3, and M 2 = 2[Ek / '(pL 2 c 2 )} introduced in §2.2. The 
well-known correlations between cloud properties ( [Larson ( 1981 )| ) can be converted into 
relationships between these parameters, too. For example, linewidth/size and density/size 
scalings of the form a v oc L 1 ^ 2 and p oc L^ 1 reported either for whole GMCs ( jSolomonj 
et al (1987)| ) or for clumps within GMCs flBalty et al (198Tf ) can be combined to show 
that cloud kinetic energy (Ek oc M 2 ) and gravitational energy (Eq oc n 2 ) scale together; 
in terms of the Mach number J\A and Jeans number rij, the observed relationship can be 
expressed quantitatively as M. = cr v /c s = l.lnj. With this observed scaling for sonic Mach 
number, the corresponding Alfven Mach number would satisfy Ai A = ct v /va = 1.7nj(3 l l 2 . 
As a consequence, these observed clouds would be magneto- Jeans unstable (nj/? 1//2 > 1) 
when M.a > 1-7, and magnetically supercritical (njP 1 ^ 2 > 1/ir) when J\4a > 0.54. We can 
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also define a fluid flow crossing time t cross = L/a v ; using the observed scalings, this implies 
t = 1 If 

Lg ±. I l cross- 

In modeling the evolution of decaying turbulence in idealized magnetized, self- 
gravitating clouds, our survey of numerical experiments brackets the range of observed 
properties wherever possible. The survey covers a variety of different initial conditions 
by varying three parameters independently: we set the size of the box to rij = 2, 3, or 
4 Jeans lengths (see eq.0), we set the ratio c 2 /v 2 A = (3 — 0.01,0.1,1, or 10 6 (see eq. 0), 
and we set the initial kinetic energy to Ex^nit/ (pc 2 L 2 ) = 25,50, or 100, corresponding to 
initial sonic Mach numbers of Ai = 7.1, 10, or 14. To allow for the decay of turbulent 
energy, we have set the range of initial kinetic energies slightly higher than estimates of 
observed cloud kinetic energies in our range of nj. The initial Alfven Mach number 
M A = (v 2 /v 2 A y/ 2 = (2(3E K ) 1/2 = Mp l/2 is 7.1, 10, or 14 for the (3 = 1 models; 2.2, 3.2, or 
4.5 for the (3 = 0.1 models; and 0.71, 1, or 1.4 for the (3 = 0.01 models. Thus, all of the 
models with (3 = 10 6 , 1,0.1 have initially significantly super-Alfvenic flow, while the models 
with (3 = 0.01 have initially sub-Alfvenic or slightly super-Alfvenic flow. For (3 = 1, all of the 
model clouds would be unstable by the magneto- Jeans criterion (if initially unperturbed), 
and also supercritical. For (3 = 0.1, the rij = 4 models would be magneto- Jeans unstable and 
supercritical, while the rij = 3, 2 models would be magneto- Jeans stable but supercritical. 
For (3 = 0.01, the nj = 4 models would be magneto- Jeans stable but supercritical, while 
the rij = 3, 2 models would be magneto- Jeans stable and subcritical. The set of parameters 
used for the simulation suite is given in the left three columns of Table [l|. 



3. Results 

3.1. Overall Cloud Energetics 

Figure 1 shows the evolution of the total "perturbed" energy E tot = Ek + SEb + Eg 
for all of the cloud models over the course of the simulations. From Figure 1, it is clear 
that low-magnetization cloud models become gravitationally-bound - defined as reaching 
Etot < - at significantly earlier times tund than high-magnetization models, and that 
the models with the same values of (3 and nj (i.e. same mean magnetic field and size) 
which become bound all do so at close to the same time, regardless of the initial level of 
turbulent kinetic energy. In Table 1, we list the time tu n d (or lower limit) at which the 
model clouds become bound. The nearly-unmagnetized (f3 = 10 6 ) cloud models all become 
bound at approximately the same time, 0.5t 9 , regardless of the initial level of kinetic energy 
or cloud size (i.e. nj). This is only slightly larger than the time (~ 0.3t g ) that a cold, 
initially- quiescent cloud would take to collapse; thus, turbulent energy does not in itself do 
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much toward preventing cloud collapse. The (3 = 1 models become bound at times between 
0.5 — l.ltg, with the larger (nj = 4) clouds collapsing before the smaller (nj = 2) ones. 
Thus, weakly-magnetized model clouds which have Alfven Mach numbers larger than one, 
and are magneto- Jeans unstable, collapse within a factor of two of the time for collapse 
of completely unmagnetized clouds. The /3 — 0.1 model clouds all have initial Alfven 
Mach numbers > 1, as well, and are all supercritical. Of these, the rij = 4 models, which 
are magneto- Jeans unstable, collapse at approximately the same time (~ O.Qt g ) as the 
nearly-unmagnetized {(3 = 10 6 ) and weakly-magnetized (3=1 models; the nj = 3 models 
(which are just barely magneto- Jeans stable) remain unbound until 0.7 — 1.0t g , while the 
magneto- Jeans stable (but supercritical) nj = 2 models evolve to 1.5t g before undergoing 
rapid gravitational collapse. For (3 = 0.01, all cloud models with nj = 2, 3 (which are 
subcritical) remain unbound up to times at least t g , with the models with higher initial Ek 
lasting until the simulations were halted at 1.5i g ; the supercritical nj = 4 clouds lose energy 
more quickly and collapse at ~ 0.8t g . 

Consistent with expectations, the gravitational energy of the subcritical cloud 
models {(3 = 0.01 and nj = 2, 3) never exceed the energy of a self-gravitating thin sheet 
Ecsheet ~ — [( n j 7I " 2 /6) — 1], over the course of the simulations (t > t g ; see Table 1). Since the 
time for an initially-quiescent, cold cloud to collapse to a sheet is just 0.25t g , it is clear that 
the time-dependent magnetic field perturbations have delayed the collapse along the mean 
field, in subcritical clouds. For the other models that do become bound, the gravitational 
energy exceeds the energy of a self-gravitating sheet at approximately the same times 
(tcondense) as the times when they become bound (hind) (see Table 1). 

For all models, Figure 1 shows that there is an initial phase of rapid kinetic energy 
loss, with the dissipation rate higher in the models with the larger initial kinetic energy 
and weaker mean magnetic fields. The time dependence of the kinetic energy dissipation 
is isolated in Figure 2, which shows the evolution of Ek for all of the cloud models on a 
log-log scale. Over the period before the matter becomes strongly self-gravitating, after a 
transient phase (during which much of the initial energy is lost) the kinetic energy declines 
approximately proportional to t -0 3 in the (3 = 0.01, 0.1 models, and proportional to t -0 6 
in the (3 = 1 models. The kinetic energy evolution in the f3 = 10 6 models does not have a 
significant period of power-law decayf] 



1 In recent work by ourselves (see Stone, Ostriker, fc Gammie (1998) ) and another group ( MacLow et al 



1998) ) , numerical simulations of decaying magnetized and unmagnetized turbulence in 3D have shown that 



the kinetic energy decline during the power-law phase is a steeper function of t. We believe that the difference 
with the present 2.5D simulations arises because motions in the symmetry direction (£) are dissipated less 
easily. Because of the difference between 2.5D and 3D, we regard the present dissipation rates as lower limits; 



- 12 - 



Since the simulations pesented in this paper all describe models of decaying, rather than 
quasi-steady-state, turbulence, there is no single number which can be used to characterize 
the turbulent decay rate for each model. However, a useful way to quantitatively 
characterize the differences in the evolution of the kinetic energy for model clouds of 
different mean magnetization and size is to compare the total kinetic energy after one 
crossing time with the initial kinetic energy. We define a characteristic crossing time for the 
flow as t cross = L/(2E K , imt ) l/2 ; for E K , inU = 25, 50, 100, t cross /t s = 0.1414, 0.1, 0.0707, 
respectively. In Table 1, we list the values of the ratio EK(t cross )/ Ek, init for each 
model cloud. For the (3 = 0.01 models, this ratio is in the range 0.42 — 0.48, with the 
value depending more on the initial kinetic energy than on nj because self-gravity does 
not become important for any of these model clouds within t cross . For the (3 — 0.1 
models, all of the models which do not become strongly self-gravitating within t cross 
have Exitcross) / Ek, init smaller than the corresponding values for the (3 = 0.01 models; 
the range is E K{t cross) / Ek, init = 0.32 — 0.46. All of the (3 = 10 6 models have smaller 
values of E K (t cross )/ E Kt i n it (in the range 0.26 — 0.32 for uncollapsed models) than the 
corresponding models with (3 = 1, and most of the (3 = 1 models have smaller values for this 
ratio (0.33 — 0.39) than the (3 = 0.1 models. Although the trends in the relative levels of 
dissipation are as expected, with the lower-/? cloud models experiencing less dissipation than 
the higher-/? models, there is less than a factor two enhancement of kinetic energy at times 
t cross in the former compared to the latter. Thus, magnetic fields do reduce dissipation, but 
when the Alfven Mach number is of order 0.5 or larger, as for the present set of models 
(and as expected from estimates of field strengths in self-gravitating molecular clouds), the 
reduction in dissipation is not large. 

In Figure 3, we show the evolution of the perturbed magnetic energy 5Eb for all 
models. Comparison with Figure 2 shows that the perturbed magnetic energy in the 
(3 = 0.01,0.1, and 1 models typically increases until its value is 1/3 — 2/3 that of the 
kinetic energy, and then both decline together. The time to reach the peak magnetic energy 
increases with increasing (3 and decreasing Ek, init, and lies in the range 0.1 — 0.25^ where 
the Alfven crossing time t& = L/v A = t s f3 x l 2 . Overall, the peak values of 5E B decrease with 
increasing (3, and depend only weakly on rij for any value of E K , init- A fit to the peak value 
of 5E B good to ±30% is 

5E B , max = 0.3E K s mit (3-°- 18 ; (9) 

this fit includes only those models with (3 = 0.01, 0.1, 1 for which 5E B reaches a maximum 
and then declines. In Table 2, we list the ratio 5E B /E K after one crossing time. For 



however, we note that the quantitative differences in the total magnitude of dissipation between 2.5D and 
3D simulations performed to date are not large (see §4.1) 
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(3 = 0.01, 0.1, and 1, this ratio varies between 0.25-0.62; for f3 = 10~ 6 this ratio is just 
4 — 10 x 10~ 6 . Because the peak value of SEb, the time at which the peak is attained, and 
the rate of kinetic energy dissipation all depend on (3, the ratio 5E B /E K (t cross ) varies with 
f3 for fixed Ek, init] the variation with nj at fixed (3 and Ek, init is more moderate. Defining 
f3 pert = cl/5v\ = 1/(2SEb), we note that for all the models with (3 = 0.01, 0.1, and 1, the 
minimum value of (3 pert is 0.1 or smaller; thus, there is dynamical growth of the perturbed 
magnetic field due to turbulence even when the initial mean field is relatively weak. The 
amount of growth of the perturbed magnetic field likely depends on the degree of helicity 
in the initial velocity field (cf. [Menneguzzi, Frisch, fc Pouquet ( 1981 )| , [Hawley, Gammie, fc| 



Balbus (1996)|) , which is close to maximal for our choice of a divergence-free initial velocity 
field. 



3.2. Evolution of Cloud Structure 

To illustrate the growth of density structure in model clouds, Figure 4 shows density 
snapshots for the model with (3 = 0.01, nj = 3, and Ek, init = 50 at time intervals 0.1t g 
(~ lMyr). This model is magneto- Jeans stable, and marginally subcritical. The early 
structure, consisting of many small clumps and interconnected filaments, is characteristic 
of the early appearance of all the models. Over time, the small structures agglomerate 
into larger structures as kinetic energy is lost to dissipation and gravity does its work. 
This general process occurs in the other models, as well. However, the late-time structural 
evolution of the models differs qualitatively depending on the values of (3 and nj. In 
subcritical models, over time gravity causes matter to slide along the magnetic fields and 
collect into larger filaments (actually the projections of sheets in this 2D geometry) which 
lie preferentially perpendicular to the mean field. Figure 4 shows how these sheets may 
continue to oscillate in time. In supercritical models, on the other hand, the self-gravity 
of the condensations that form is sufficient to overwhelm magnetic stresses, leading to the 
growth of rounder clumps and eventual gravitational collapse. 

Figure 5 compares the development of structure in three cloud models which differ 
only in the magnitude of their mean magnetic field - having (3 = 0.01,0.1, and 1. The 
(3 = 0.01 model is magneto- Jeans stable and (marginally) subcritical; the (3 = 0.1 model 
is magneto- Jeans stable but supercritical; the (3 = 1 model is magneto- Jeans unstable 
and supercritical. The model clouds all have the same size nj = 3, and all have the same 
value of Ek, init = 50 (and, in fact, identical initial spatial velocity perturbations). In the 
(3 = 0.01 model, the magnetic field lines remain close to their original, parallel configuration; 
early-time large oscillations of the matter clumps and filaments bend and compress the field 
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locally. At the final time shown (t = 0.8t g ), the matter density remains less than p/p— 100 
everywhere. In the /? = 0.1 model, the magnetic field is less rigid and hence there are much 
larger initial departures from parallel lines. Local regions with density log(p/p) > 2 appear 
in the t/t g = 0.5 snapshot (these contain 3.6% of the total mass), and grow in volume and 
mass so that at t/t g = 0.8, 18% of the mass is at log(p/p) > 2. The magnetic field becomes 
more uniform as kinetic energy is lost, but there remain large departures from straight 
lines in the high-density regions. In the (3 = 1 simulation, the magnetic fields immediately 
become strongly bent and stretched, and by t/t g = 0.8 has even undergone noticable 
reconnection in the high-density region (due to numerical effects). The matter becomes 
strongly clumped, with 8 % of the mass at log(p/p) > 2 at t/t g = 0.5, and 62 %(55%) 
of the mass at log(p/p) > 2(3) at t/t g = 0.8. In this last frame, a single zone contains 
47% of the total mass; because the density in this zone is high enough to make the Jeans 
length smaller than the grid scale, the simulation is not well-resolved at this point. The 
(3 = 10 6 simulations (not shown) are also initially filamentary, and gradually coagulate due 
to gravity. The late-time structures consist of rotationally-supported condensations (with 
trailing spiral arms); this sort of structure does not develop in the lower-/? models because 
magnetic stresses redistribute angular momentum away from gravitationally-collapsing 
regions. 



As described in §3.2, the overall structure in our model clouds evolves in time due to 
the combined action of gravity, Reynolds stresses, and magnetic stresses. At the same time 
as the overall structure evolves, the detailed distributions of matter in the models evolves as 
well. These distributions can be characterised in a variety of ways. We find it informative to 
consider two related distributions: the fraction of volume dV/Vtot, and the fraction of mass 
dM/M tot = pdV I (pV tot ) , as a function of the logarithmic density enhancement /decrement 
log(p/p). In Figure 6, we show how these volume and mass distributions appear for the 
cloud model snapshots portrayed in Figure 5. The values of the volume-averaged and 
mass-averaged density logarithms 
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Statistics of Matter Distributions 



(log(p/p))y 
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(log(p/p)) M 
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\og(p/p)pdV(p) 
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(11) 



are also indicated, for each distribution. 



- 15 - 



A number of common features are apparent for all the models shown in Figure 6; these 
trends apply to our other models, as well: (a) The mass distributions are offset toward larger 
density than the volume distributions; this simply means that matter is clumped, with 
the mass-averaged density logarithm larger than the volume-averaged density logarithm. 
Because most of the volume is at lower-than-average density, while most of the mass is at 
higher-than-average density, (log(p/p))y is negative and (log(p/p)) M is positive, (b) At 
earlier times, the distributions are very roughly log-normal in shape, while large tails evolve 
in time, (c) The values of (log(p/ p))v and (log(p/p))M are nearly equal in magnitude, 
at early times. It is easy to show that for a lognormal distribution, this will be true by 
definition, (d) Over time, the mass and volume distributions become more separated, for 
any given model. This is due to the action of gravity increasing the clumping of matter, (e) 
At the earliest time, before self-gravity becomes important, the /3 = 0.01 model has greater 
density contrast than the j3 — 0.1, 1 models. 

Given the strong temporal evolution in any of our models, no single distribution 
characterises its density structure. Because the density distribution in a real cloud is a 
fundamental property, however, it is important to investigate in our model clouds how the 
density distribution depends on the input parameters. If a clear correlation exists between 
the density distribution and input parameters in our model clouds, then the density contrast 
in real clouds (estimated from molecular mapping), and the more detailed distributions of 
column density (determined from extinction mapping), can potentially be used to constrain 
other cloud properties such as the mean magnetic field strength. While the limitations 
of our isothermal equation of state may affect our conclusions in detail, we still consider 
it important to investigate how turbulence and magnetic fields alone impact the density 
structure that develops. 

To isolate the effect of the magnetic field in influencing the density structure in our 
model clouds, we have chosen to analyze the models early in their evolution, before the 
effects of gravity become significant. Since the velocity dispersion - corresponding to the 
kinetic energy in our models - is the primary direct observable in real clouds, we select 
instants in each model cloud's evolution when the kinetic energy has a fixed value; because 
Ek decays at different rates for different /3, these analysis times differ for the different cloud 
models. To make quantitative comparisons among the models, we evaluate the volume- 
and mass- averaged density logarithms at times in each model when the kinetic energy has 
decayed to half of its initial value. Thus, for models with initial E K:init = 100, 50, 25, the 
respective values of the Mach number at the time of density structure measurement are 
A4 = 10,7.1,5. The computed values for (log(p/p))y and (log(p/p))M are listed in the last 
two columns of Table 1. In a schematic picture of clouds as consisting of denser clumps 
embedded within a more tenuous interclump medium, we may think of (log(p/p))M as 
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describing the magnitude of the positive density contrast in the "clump" component, and 
(log(p/p))y as describing the magnitude of the negative density contrast in the "interclump 
medium" . 

Based on the computed values of (log(p/p))y and (log(p/p))M tabulated, we can make 
a number of observations. These remarks apply to the set of 27 models with (3 = 0.01, 0.1, 1; 
since the unmagnetized models have unrealistically low field strength, and generally become 
self-gravitating earlier, we do not consider them here. At fixed Ex and (3, we find the 
variation of the the density contrast with nj is (except for the (3 = 0.01, Ek, % n it — 0.25 
models) less than 0.15; thus, these variations are not principally due to gravitational effects 
(as desired). We also find that for fixed (3, the density contrasts (log(p/p))y,M generally 
increase logarithmically with increasing E K . Finally, we find that at fixed E K , all except for 
two models with (3 = 0.01 have larger density contrasts than the corresponding (3 = 0.1, 1 
models at all raj. In general, we find that the variation of density contrast with Mach 
number at fixed f3 is smaller than the variation of density contrast with f3 at fixed Mach 
number, for the parameter range studied. 

In Figure 7, we plot the values of the clump and interclump density contrasts 
vs the total magnetic energy density, so that the models with different mean (3 fall 
in three separate groups. Within each group, increasing (perturbed) magnetic energy 
correlates with increasing kinetic energy, accounting for the increase shown in the 
density contrast toward larger total Eb- The (3 — 0.1 points lie quite close to the lines 
|(log(p/p))v iA f| = C\og(B 2 /(8npc 2 s )) with C = 0.2 - 0.25, and the (3 = 0.01 points lie near 
the extension of these lines. The (3 = 1 points all lie well above these lines, and typically 
have density contrasts that fall between the (3 = 1 and (3 = 0.01 models, at a given E K . 

Physically, the decrease in the density contrast in our models in going from (3 = 1 
to (3 — 0.1 at fixed Ek can probably be attributed to a increase in the effective pressure 
provided by magnetic fields: for a fixed value of the Mach number = (2EkY^ 2 , the 
Alfven Mach number decreases proportional to /5 1//2 , and ram-pressure induced density 
contrast should decrease as the mean magnetic field strength increases. The increase of 
the density contrast in going from (3 — 0.1 to (3 — 0.01, on the other hand, is likely the 
signature of the importance of transverse magnetic fields in compressing the medium. At 
quasiequilibrium kinks in a strong magnetic field, matter is compressed until the central 
thermal pressure pc 2 equals the external magnetic pressure Bj_/8n. For a nonlinear 
transverse wave, Bj_/4ir ~ pv\, implying p/p ~ (l/2)(fj_/c s ) 2 at kinks (independent of 
(3). In this situation, the induced density contrast would not decrease with increasing 
mean magnetic field strength. The measured tendency for matter to collect more strongly 
at nulls of the transverse magnetic field for lower-/? simulations may in part owe to the 
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increasing virulence at low (3 of the decay instability (in which Alfven waves pump the 
growth of field-aligned compressive waves and backscattered Alfven waves; e.g. |Sagdeev""fc" 



Galeev (1969)| , |Goldstein (1978)| ), which may cause the lowest-/? models to approach the 



kinked-field state more rapidly. 



Our results for density contrast can be compared with those reported by Padoan, 



Jones, fc Nordlund ( 199771 , an d those of [Passot fc Vazquez- Semadeni (199"8")| for isothermal 



models. The former cites (based on 3D simulations of [Nordlund fc Padoan (1998]| ) a result 



-0.51og(l ±.M 2 /4) for the value of (log(p/p))y, whereas Passot fc Vazquez- Semadem 



1998)| (from ID unmagnetized simulations) find the result — 0.22.M 2 for the same quantity. 



The measurements in our Table 1 are made for A4 = 10, 7.1, 5, so that the padoan, Jones 



fc Nordlund (1997)] formula predicts ±0.7, 0.6, 0.4, while the |Passot fc Vazquez-Semadcn" 



(1998)| formula predicts ±22, 11, 5 for the values of (log(p/ p))m,v- Because |Padoan, Jones 



Nordlund ( 199*7] do not state the value of (3 (or other properties of their simulations), we 



cannot make direct comparisons with our results; nevertheless, the range of values ±0.2 — 1 
for the density contrast in Table 1 is consistent with the values they predict. Our range 
of values is considerably smaller than the values predicted by |Passot fc Vazquez-Semadeni 



1998)| for purely ID turbulence. We have also recently found, using forced turbulence 
simulations in 3D ( |Stone, Ostriker, fc Gammie (199"8"j| ), comparable values for the density 
contrast to those obtained in the present 2 1/2 D models. 

As a cautionary remark, we note that our conclusions may potentially be affected 
by the assumption of isothermality. Passot fc Vazquez-Semadeni (1998]] found, using ID 
simulations, that a power law may develop over a range of densities larger (smaller) than 
the mean for polytropic 7 < 1 (7 > 1); although they do not state how (\og(p/ p)) m,v 
depend on 7, presumably these quantities vary as well. The power-law regime increases 
with increasing M. and |1 — 7). However, because the ID simulations at 7 = 1 appear to 
overestimate the dispersion and mean density contrast relative to 2.5D and 3D simulations 
at any value of M. and 7 the range of values of log(p/p) over which a power law would 
occur in higher-dimension simulations is likely much smaller than those presented by Passot 



fc Vazquez-Semadeni (19981 . 111 addition, since 7 ^ 1 (7 < 1) probably applies for densities 



larger (smaller) than ~ten times the mean for molecular cloud conditions (cf. |5calo et~al 



1998)|) , it is not clear whether a power law would develop within the molecular component 



for a substantial portion of a real cloud's internal density distribution. f\ At the time when 
density statistics are computed in our models, most of the material lies at < log(p/p) < 1, 



5 On the other hand, in the large-scale ISM for which the effective 7 < 1 over a larger range of densities 
across different gas phases, a larger power-law regime should develop; this result has been found in the 
analysis by Scalo et al (1998) of non-polytropic 2D simulations. 
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so that for any realistic 7 the temperature would not have departed far from the mean 
value (as is indeed observed in molecular clouds); in this sense, our isothermal assumption 
is self-consistent for the simple mean density contrast statistic which we have measured. 
Thus, while calculations incorporating heating and cooling are clearly warranted, we 
speculate that at least the values of (log(p/p))jv/ which would be obtained would not be far 
different from those calculated here with an isothermal assumption. We discuss potential 
implications of our results on density contrast in §4.3. 



4. Discussion 

4.1. Summary and comparison with other work 

In this paper, we present the results of a suite of 2.5-dimensional simulations following 
the evolution of idealized cloud models with dynamical parameters comparable to those 
believed to hold in Galactic molecular clouds. Our models, which assume cold, isothermal 
conditions, are initiated with random, turbulent velocity fields, and include the effects of 
self-gravity and magnetic forces. We parameterize the importance of magnetic fields by 
(3 = c^/v\ (see eq. [J), and the importance of self-gravity by the Jeans number rij = L/Lj 
(see eq. ||). We follow the evolution of cloud models over time as they dissipate kinetic 
energy in shocks, and as the clumps, sheets, and filaments formed by the action of magnetic 
and Reynolds stresses coalesce under their own self-gravity. 

In §3.1, we compare for different models the times needed for clouds to become 
gravitationally bound, and the times needed to reach a state with stronger gravitational 
binding than a thin sheet. We find that only the models which are magnetically subcritical 
((3 l l 2 nj < 1/ir) survive gravitational collapse for times significantly longer than t = t g \ 
in §4.2 below we comment on the implications of these results. We have also compared 
the evolution of kinetic energies, and perturbed magnetic energies, among our models. In 
examining the amplification of the perturbed magnetic field, we found that the peak value 
of SEb attains higher amplitudes, and occurs at earlier times, for the lower-/? models (see 
eq. [|); at the peak, the value of SE B is 1/3 — 2/3 times the value of E K . To quantify 
the decay of turbulent energy, we compare the kinetic energy after a cloud flow crossing 
time (i.e. an "eddy turnover time," in the terminology of incompressible turbulence) 
with the initial kinetic energy. We find, for models that have not become strongly 
self-gravitating at this time, that the kinetic energy is a fraction 0.26-0.29, 0.33-0.36, 
0.32-0.46, and 0.42-0.48 of the initial value, for (3 = 10 6 , 1, 0.1, and 0.01, respectively, so that 
(Ek) /(AEk/ At) = 0.9 — lAt cross . Thus, while the models with stronger fields have lowered 
dissipation, the effect occurs at a level less than a factor of two, for realistic field strengths. 
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Although the 3D simulations of [MacLow et al (1998) and Padoan fc Nordlund (1997) have 
somewhat different parameters, they also find kinetic energy reduced to 0.3-0.5 of the 
initial value, after one crossing time, and differences of no more than a factor two between 
their low-/? and high-/? cases. In recent forced, 3D simulations ( [Stone, Ostriker, fc Gammic 
1998)|) , we have found dissipation times E\y/E for the total perturbed energy a factor 



1-1.75 times smaller than the loss time for half of the wave energy for the present 2 1/2 D 
models; for purely 2D models, the dissipation rate is even higher than in 3D. On the other 
hand, the dissipation in our set of ID models (Paper I) was much less than in the present 
set, probably because the ID simulations' initial velocity fields were purely transverse, and 
the model symmetry inhibited dissipation. Since real clouds are more similar to the 2.5D 
and 3D models than the ID models, we conclude that realistic turbulent dissipation times 
E K /E K are not likely to be more than twice the cloud fluid crossing times - i.e., £ 10 
Myr - for any likely value of the mean magnetic field. Although the hope of substantially 
lowering turbulent dissipation rates in molecular clouds by including magnetic fields ( [Aronsj 

Max (1975)| ) may not be fully realized, the feedback provided by star formation (e.g. 
|Norman fc Silk ( 1980 )| ) on timescales comparable to the turbulent decay (see §4.2) may 
yet allow a quasi-steady state to exist (e.g. [McKee (1989)| , Passot, Vazquez-Semadeni, fc 
Pouquet (19951 , IGammie fc Ostriker ( 19961 ). 

We present, in §3.2, snapshots showing density contours, magnetic field lines, and 
the velocity field in cloud models with varying levels of the mean magnetic field strength 
(i.e. f3). All models show highly inhomogeneous density structure well before gravity 
becomes important in the dynamics. However, as gravity becomes important, the matter 
in subcritical clouds collapses into sheets, whereas in supercritical clouds it collapses into 
rounder clumps (and reaches much higher densities). In §4.2 below, we summarize our 
results (based on the analysis in §3.3) on the distributions of matter in our models, and 
comment on implications for interpreting the clumpy structure of observed clouds. 



4.2. Implications for cloud support and the initiation of star formation 

Observed molecular clouds have column densities in the range Nh 2 = 3 — 10 x 10 21 cm -2 
(e.g. [Larson ( 1981 )| , |Myers fc Goodman (1988"b~J| , [Blitz (1993)| ). The corresponding magnetic 
field strengths required for such clouds to be cross-field ( "magneto- Jeans" ) stable (see §2.2) 
are in the range 7 — 24/iG, within the range of field strengths that are believed to prevail 
on the large scale in molecular clouds (§1). However, typical observed field strengths fall 
short of the value 23 — 77/iG which would be required to render the same whole clouds 
magnetically subcritical; i.e., able to resist gravitational fragmentation if all of the cloud's 
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mass is collected into a sheet, rather than homogeneously filling a volume with sides 
comparable to its lateral extent. We note that projection effects tend to decrease the 
observed (line-of-sight) field strength relative to the 3D value, and increase the column 
density relative to its value along the shortest axis, so that observed B, Nh 2 values might 
still be consistent with marginally critical clouds. 

In the simulations presented in this paper, we show that regardless of the initial 
turbulent energy content of our model clouds, only those clouds which are magnetically 
subcritical remain uncollapsed throughout the simulations for times exceeding 1.5t g - 
corresponding to m 15Myr for typical conditions. In fact, the only supercritical clouds 
that survive uncollapsed beyond £ g ~ lOMyr are relatively small (nj = 2, L ~ 4pc) and 
magneto- Jeans stable. Very low magnetization clouds collapse in 0.5t g ~ 5Myr; this puts a 
lower limit on the collapse time for supercritical clouds. If the Milky Way molecular clouds 
containing most of the mass have comparable field strengths to those in observed regions, 
and column densities are at the high end of the range mentioned above, then they would be 
supercritical and our simulations would imply that they cannot survive for periods longer 
than 5 — lOMyr without some of part of their interiors becoming strongly gravitationally 
bound, and, we expect, beginning to form stars. Alternatively, if most of the molecular 
material is in low column density, subcritical clouds, then the dissipation of turbulence 
incorporated from formation would lead to a flattened configuration supported by magnetic 
fields (if there is no turbulent resupply from outside sources), but not dependent on 
turbulent input by stars to avoid further collapse. This is one of the main conclusions of 
the present work. 

Observationally-based estimates of the time to initiate star formation are necessarily 
indirect, and subject to major uncertainty. Based on |Williams fc McKee (1997)| 's analysis 



of the joint distributions of OB associations and GMCs (assuming all molecular material 
is associated with discrete clouds with a spectrum of masses dN/dM oc M -L6 ), if star 
formation is equally distributed among all clouds, then roughly 80% of the Galactic 
molecular material would reside in clouds containing one or more OB stars. If we take 
this fraction as representing the portion of a cloud's life during which it contains an OB 
association, and use 20 Myr (piaauw ( 1991 )|) as the typical observed lifetime of associations, 
this implies that associations form in clouds after 20% of their 25 Myr lifetimes, essentially 
as soon as our models predict that star formation is initiated in supercritical clouds (5-10 
Myr). On the other hand, even in the largest existing cloud catalog ( [Solomon et al (1987)| ) 
created for the inner Galaxy, ~ 60% of molecular material was not actually assigned to 
discrete clouds, perhaps because it has not been sufficiently warmed internally by massive 
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stars to stand out above the background ( [Solomon fc Rivolo ( 1989)1) .p| If this cold material 
is in clouds with the same mass distribution as the warm clouds, we would infer that 60% 
of a cloud's lifetime passes before OB stars form (and it becomes warm). Under this second 
scenario of unequally-distributed OB stars, typical cloud lifetimes would be 50 Myr, and 
a period of 20 — 25Myr must therefore be spent in forming exclusively low mass stars, if 
the clouds are supercritical and the first stars form after 5 — lOMyr. If, on the other hand, 
much of the cold molecular material is distributed in small, subcritical clouds that never 
make OB associations, then OB star formation in supercritical clouds (comprising less than 
half of the total molecular material) would again be inferred to commence soon after clouds 
first begin to have portions of their interior collapse. 

Based on the above discussion, the range of star formation scenarios permitted by 
observations is fairly large at present. In particular, the crucial question of the fraction 
of molecular material residing in sub- and super-critical clouds is essentially unknown. 
However, the planned high-resolution Boston University/FCRAO Galactic Ring Survey 
in 13 CO and CS (M. Heyer, personal communication) should considerably constrain 
the molecular cloud distribution function, and correspondingly constrain the possible 
evolutionary scenarios for star formation and cloud destruction. This, in turn, should allow 
us better to evaluate the conclusions of this paper on the time to initiate star formation in 
magnetized GMCs. 



4.3. Implications for clumpy structure of GMCs 

The other finding which we would like to highlight concerns the distribution of density 
found in our simulated clouds, and the potential to use observed density distributions in 
real clouds as a diagnostic of the magnetic field strength. We find, in concurrence with 
other recent work on isothermal cloud models, that the shapes of the density distribution 
functions before self-gravity becomes important are roughly log-normal. Our own analysis 
of density distributions presented in §3.3, however, primarily focuses not on the shape 
of the density distribution function, but instead on the values of volume-weighted and 
mass-weighted means of the logarithmic density, (log(p/p))y and (\og(p / p)) m ■ These 
represent the magnitudes, respectively, of the (negative) density contrast in the tenuous 
interclump region, and the (positive) density in the cloud clumps, relative to the mean 



6 For a fixed number of OB associations, reducing the fraction of clouds containing massive stars increases 
the photoevaporation rate for those clouds that do contain associations, making estimates of cloud destruction 
timescales more consistent with observed OB association lifetimes ( Williams fc McKee (1997) ). 
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cloud density. For a lognormal distribution, (log(p/p))y and (log(p/p))M will be equal in 
magnitude and opposite in sign; this result holds for most of our models, when analyzed 
at early times. We find that the values of the mean density contrast depend both on 
the amplitude of the turbulence (i.e. the value of Ek), and on the strength of the mean 
magnetic field and its perturbations (i.e. the value of Eb). Specifically, we find that 
for a given j3, the density contrast increases with Ek, and that virtually all of the more 
strongly-magnetized (/3 = 0.01) models have larger mean density contrast than the more 
weakly-magnetized {(3 = 1,0.1) models with the same value of Ek- We also find that the 
(3 = 1 models generally have higher density contrast than the /? = 0.1 models. Given the 
limitations of the current set of simulations, we did not attempt to formulate a more general 
theory for the relationship between the mean density contrast and the values of Ek and 
(3. We defer this much-desired goal until the completion of surveys of 2D and 3D forced, 
non-self-gravitating MHD turbulence which are now underway. With 3D simulations, one 
can compute column density distributions, for better comparsion with directly-observed 
cloud properties. We emphasize that our model clouds are highly idealized, and in particular 
we employ an isothermal equation of state which may affect our conclusions (see §3.3.). 
Still, it is significant theoretically, and potentially useful astronomically, that the mean 
magnetic field strength affects the density distribution. 

While follow-up studies are clearly needed, in particular considering non-isothermal 
models, it is nevertheless tempting to compare our present results with observations. 
Molecular clouds are often characterized as consisting of clumps embedded in an interclump 
medium, with a fraction 9h = Mn/Mtot of the mass, and fu = Vn/Vtot of the volume 
contained in the high density component (ph/p = &h/ fa), and the remainder of the mass 
and volume dispersed in a low-density component (pl/p — (1 — @h)/(1 — fu))- From 
13 CO maps of Orion A and the Rosette clouds, Pally et al (1987)] and [Williams, Blitz, ~Ez\ 



|Stark (1995)| estimate 6 H = 0.75,0.77 and f H = 0.1,0.08, respectively. With this simple 
clump/interclump structure, we use 

(\og(p/p)) v = f H \og{ PH /p) + (1 - f H ) hg{p L /p) (12) 

and 

(log(p/p)) M = 6 H \og(p H /p) + (1 - H ) \og(p L /p). (13) 

to compute density contrast magnitudes of (log(p/p))y = —0.4,-0.5 and 
(log(p/p))M = 0.5,0.6 for the two clouds. The clumps defined in these observations 
are not believed to be self-gravitating; thus, it is fair to compare to the model density 
contrast magnitudes computed in §3.3. Of the present models with (3 = 1,0.1,0.01 
(corresponding to B ~ 1.4,4.4, 14.4/iG), Figure 7 shows that only the (3 = 0.01 simulations 
with Mach number 10 have density contrasts as large as ±0.5; thus, the stronger-field 
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cases seem to be in better agreement with observations, when taken at face value. As 
argued in §3.3, both (3 > 1 and /3 < 0.01 could yield larger density contrast, although 
the former (corresponding to sub-/zG fields) is unlikely to exist in the present-day Galaxy. 
The implication that Bq £ 15/xG is far from definitive, but demonstrates the potential 
to develop magnetic field diagnostics from other cloud properties such as the inferred 
clump/interclump density contrast (which is albeit somewhat model-dependent), and the 
observed distributions of column density (obtained from extinction studies, and requiring 
only an Ay — > Nh 2 conversion). We view this as an an important direction for future 
research. 

We are grateful to B. Elmegreen, J. Ostriker, F. Shu, and the referee J. Scalo for careful 
reading and helpful comments on the manuscript of this paper, and D. Neufeld for useful 
discussions on molecular cooling. 
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Table 1. Evolutionary Characteristics of Model Clouds 
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h t C ondense is defined as the time at which E G /E GyShcct = 1 where E Gy shcc t = -[(n 2 T 7r 2 /6) - 1]. 
c Evaluated at the crossing time t cr = L/(2Ek, init) 1 ^ 2 

d Volume-weighted average of \og(p/p), evaluated at the time when Ek — EK,init/2 (see text). 
c Mass- weighted average of log(p/p), evaluated at the time when Ek — EK,init/2 (see text), 
t Extrapolated 



* Self-gravity has become unresolved at the grid scale before this measurement is made 
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Fig. 1. — Evolution of the total specific energy E tot /cf = [Ek + SEb + Eg) / for all models. 
Solid, dotted, and dashed curves have rij = 2, 3, and 4, respectively. 




Fig. 2. — Evolution of the specific kinetic energy for all models. Solid, dotted, and dashed 
curves have rij = 2, 3, and 4, respectively. 




Fig. 3. — Evolution of the perturbed magnetic energy (per unit mass) for all models. Solid, 
dotted, and dashed curves have rij = 2, 3, and 4, respectively. 




Fig. 4. — Density snapshots of the (3 = 0.01, nj = 3, Ek, mit = 50 model cloud at temporal 
intervals of 0.1t g . Density contours are at equal intervals in log(p/p) = 0, 0.5, 1., .... 



Fig. 5. — Comparative evolution of cloud models with different magnetization. Left, center, 
and right columns have f3 = 0.01, 0.1, and 1, respectively. Ek, m%t — 50 and rij = 3 for all 
three runs. Snapshots at times t/t g = 0.2, 0.5, and0.8 show density contours (log(p/p) = 
0-0.9:dark blue, l-1.9:light blue, 2-2.9:yellow, 3-5:magenta), velocity vectors (red), and 
magnetic field lines (green). 
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Fig. 6. — Comparative evolution of density distributions in cloud models with different 
magnetization. Same parameters and intervals as in Fig. 5. The solid histograms show the 
distribution of volume dV/V to t, and the dashed histograms show the distribution of mass 
dM/M tot , as a function of density, for each snapshot. The numbers in the upper-right of 
each panel indicate the values of the means (log(p/p))M (top) and (log(p/p))y (bottom) for 
each distribution. 
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Fig. 7. — Mass-averaged (upper points) and volume-averaged (lower points) density contrast 
in model clouds when E K = 0.5E Kt init , as a function of total magnetic energy. Vertical 
dashed lines separate results from models with (3 = 1,0.1,0.01. Point type denotes Jeans 
number ny. triangles, squares, pentagons for nj = 2,3,4, respectively. The dotted lines 
indicate logarithmic slopes of ±0.2 and ±0.25. 



